Ab-inito Linear Scaling Response Theory: 
Electric Polarizability by Perturbed Projection 
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A linear scaling method for calculation of the static ab inito response within self-consistent field 
theory is developed and applied to calculation of the static electric polarizability. The method is 
based on density matrix perturbation theory [Niklasson and Challacombe, PRL (LJ9425)], obtain- 
ing response functions directly via a perturbative approach to spectral projection. The accuracy 
and efficiency of the linear scaling method is demonstrated for a series of three-dimensional water 
clusters at the RHF/6-31G** level of theory. Locality of the response under a global electric field 
perturbation is numerically demonstrated by approximate exponential decay of derivative density 
matrix elements. 

PACS numbers: 31.15.Ar,31.15.Md,31.15.Ne,33.15.Kr,36.40.Cg 
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Linear scaling methods that reduce the computational 
complexity of electronic structure calculations to O(N), 
where N is system size, impact disciplines that demand 
quantum simulation of increasingly large and complex 
systems 0, 0, & 0, IE §- To date, the most success- 
ful and prolific applications of linear scaling technolo- 
gies are ground state studies involving empirical model 
Hamiltonians. To a lesser degree, ground state applica- 
tions using ab initio models with A-scaling contribute to 
a variety of fields, typically at the Self-Consistent-Field 
(SCF) level of theory and requiring in some cases parallel 
implementations to reach levels of applicability compara- 
ble with empirical models. These SCF theories include 
the Hartree-Fock (HF), Kohn-Sham Density Functional 
(DF) and hybrid HF/DF models. 

Beyond ground state methods, little attention has been 
given to linear scaling algorithms for the computation of 
dynamic and static response properties; the latter includ- 
ing the nuclear magnetic shielding tensor Q, the rota- 
tional g-tensor indirect spin-spin coupling constant 
[3, 0| , third order properties such as the first hyperpo- 
larizability [TfT | and polarizability derivatives such as the 
Raman intensity fl2 . Il3| . Dynamic properties may be 
computed using linear scaling algorithms to propagate 
the density matrix jQ, llij in the time domain, followed 
by convolution to obtain the spectral response. In this 
way, Yam, Yokojima and Chen [15| have recently demon- 
strated linear scaling computation of the absorption spec- 
tra for one-dimensional polymers at the local density level 
of theory, but requiring ~ 14, 000 time steps. In the static 
zero frequency limit, solving the coupled-perturbed self- 
consistent-field (CPSCF) equations using standard algo- 
rithms is likewise difficult for large systems. Several al- 
gorithms have therefore been proposed for solving the 
CPSCF equations that may be capable of achieving a 
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reduced scaling better than 0(N 3 ) [3, Ell- 
in this letter, we develop the density matrix pertur- 
bation theory of Niklasson and Challacombe [18| for N- 
scaling solution of the ab initio CPSCF equations, and 
demonstrate the early onset of linear scaling for the ac- 
curate calculation of the first electric polarizability of 
three-dimensional systems using large basis sets. This 
perturbed projection method is general and can be ex- 
tended to a variety of static properties. 

Algorithms for linear scaling self-consistent-field (SCF) 
theory exploit the quantum locality of non-metallic sys- 
tems, manifested in the approximate exponential decay of 
the density matrix with atom-atom separation, through 
the effective use of sparse matrix methods and iterative 
approaches to spectral projection 0,H3|. This quantum 
locality should, in principle, extend also to the deriva- 
tive density matrices central to the CPSCF equations. 
Indeed, exponential decay of the derivative density ma- 
trix within ab initio SCF theory has been demonstrated 
numerically for local nuclear displacement However, 
standard approaches to the CPSCF equations [13)113,13 
do not admit exploitation of this locality, as they are 
based on perturbation of the wave function, requiring 
0(N 3 ) eigensolution and typically 0(N 5 ) transformation 
of two-electron integrals into the eigenbasis. Avoiding 
both eigensolution and integral transformation, Ochsen- 
feld and Head-Gordon [U and later Larsen et al. 
proposed iterative solutions to the CPSCF equations in- 
volving purely non-orthogonal representations. In both 
of these approaches, a linear system of equations contain- 
ing commutation relations must be solved. 

Recently, a formulation of density matrix perturbation 
theory has been proposed by Niklasson and Challacombe 
(NC) [13 that presents a new opportunity for solving 
the CPSCF equations within the context of linear scaling 
spectral projection [l^. l2fj| . The new approach is based 
on the relationship between the density matrix T> and the 
effective Hamiltonian or Fockian T through the spectral 
projector (Heaviside step function) V = 0(p,I — J 7 ), where 



2 



the chemical potential ft determines the occupied states 
via Aufbau filling. Spect ral proj ection can be carried out 

in a number of ways [TliflliilillHIHIilliJ, Wlth 
perhaps the best know alg orithm being McWeeny's cu- 
bic purification scheme 26] . More recently new recursive 
polynomial expansions of the projector have emer ged , 
such as the second order trace correcting (TC2) [l9j 



and fourth order trace resetting (TRS4) [2Gj purification. 
These new methods (TC2 and TRS4) have convergence 
properties that depend only weakly on the band gap, do 
not require knowledge of the chemical potential and per- 
form well for all occupation to state ratios. In the NC 
approach, the perturbation expansion is developed within 
the reference groundstate projector allowing order by or- 
der collection of terms at each iteration, establishing a 
quadratically convergent sequence for the response func- 
tions. 

We now proceed with development of the perturbed 
projection approach, which will be outlined for the con- 
crete case of perturbation by an electric field within the 
Hartree-Fock (HF) theory. However, the perturbed pro- 
jection method is general and extensible to DF and hy- 
brid HF/DF models, and to other static perturbations 
through high order. 

In the following, the indexes a, b, . . . refer to perturba- 
tion order, while . . . mark the iteration count. The 
symbols T>,J-,... are matrices in an orthogonal repre- 
sentation, while D,F, . . . are the corresponding matrices 
in a non-orthogonal basis. The transformation between 
orthogonal and non-orthogonal representations is carried 
out in O(N) using congruence transformations [32I 
provided by the AINV algorithm for co mputing sparse 
approximate inverse Cholesky factors yj, |35|, |3g . 

Within HF theory, the total electronic energy E tQt of 
a molecule in a static electric field £ is 

E tot (£)=Tr[D(h°+fxE)] + ^Tr[D(J(D)+K (£>))], (1) 

where D is the density matrix in the electric field, h° 
is the core Hamiltonian, fj, is the dipole moment matrix, 
J(-D) is the Coulomb matrix and K(D) the exact HF 
exchange matrix. The total energy may be developed in 
the perturbation expansion 



E tot (£) =E tot (0) + Y» a £ a - 
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where oi a b is the first order polarizability, \x a is the dipole 
moment and £ is the electric field in direction a. The 
polarizability is the second order response of the total 
energy with respect to variation in the electric field [2^ 
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The first order density matrix derivative T> a in the a di- 
rection is obtained by variation of both the spectral pro- 
jector 6 and the Fockian T = T a + J2 a ^ OJ70 + ■ ■ • as 
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6(fil- (T° + £ a T a )) 
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The HF Fockian F° in the non-orthogonal basis is 
F° = h° + J{D°) + K(D°), where the Coulomb ma- 
trix J may be computed in O(NlgN) with a quantum 
chemical tree code (QCTC) ,37] and the exchange ma- 
trix K computed in O(N) with the ONX algorithm that 
exploits quantum locality of D° [38j . Likewise the deriva- 
tive Fockian, F a = [i a + J(D a ) + K(D a ), may be com- 
puted with the same algorithms in linear scaling time if 
D a manifests decay properties similar to D°. A similar 
equation holds for the derivative Fockian within DF and 
hybrid HF/DF theories with addition of the exchange- 
correlation matrix V£ C (D°, D a ) 

In our approach to linear scaling computation of the 
polarizability a a b, the ground state density matrix T>° is 
computed using a spectral projection algorithm such as 
TC2 Il9l in conjunction with sparse atom-blocked linear 
algebra|2ol lioj . Linear scaling is achieved for insulating 
systems through the dropping (filtering) of atom-atom 
blocks with Frobenious norm below a numerical threshold 
(~ l(T 4 -lCr 6 ). At SCF convergence the TC2 algorithm 
generates a polynomial sequence defining the groundstate 
projector, from which expansion of the derivative density 
matrix can be obtained term by term. 

The derivative density matrix and derivative Fock- 
ian depend on each other implicitly and must be 
solved for self-consistently via the coupled-perturbed self- 
consistent-field (CPSCF) equations. The necessary and 
sufficient criteria for convergence of the CPSCF equa- 
tions are [T a , V°] + [J 70 , V a ] = and V a = V a V° + V°V a 
|4lj . Solution of the CPSCF equations with perturbed 
projection involves the steps 



F° = ^ + J{D a n ) + K(D a n ) 
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V a n+l = —e{ixI-{F° + £"?«)) 
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(5b) 

(5c) 
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with starting point Dq = 0. In step J5a|i. F£ is con- 
structed in 0{N) using the QCTC ^ and ONX [H 
algorithms in MondoSCF 01- Next, Weber and Daul's 
DDIIS algorithm for convergence acceleration of the CP- 
SCF equation [43 is used to optimize the c& coefficients in 
step l)5b[l. keeping the last m steps in the extrapolation. 
Then, the density matrix derivative is obtained in 

step (IHcll as T>n + \ = lim^oo Xf via the NC density ma- 
trix perturbation theory, based on the TC2 projector: 
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Tr[X?] > N e (6) 
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FIG. 1: Total CPU time of the fifth CPSCF iteration for 
the water cluster sequence with the 6-31G and 6-31G** basis 
sets and the GOOD and TIGHT numerical thresholds (see text) 
controlling numerical precision of the result. The lines are fits 
to the last three and four points, respectively. 
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TABLE I: Average polarizabilities a = (a xx + a yy + 
ot Z z)/3NH 2 o in a.u. for a sequence of water clusters at the 
RHF/6-31G and RHF/6-31G** levels of theory. A compar- 
ison is made between results obtained with the GAMESS 
quantum chemistry program |45| and those calculated with 
MondoSCF using different numerical approximations (see 
text) controlling precision of the linear scaling algorithms. 
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TIGHT 


GOOD 


TIGHT 


10 


4.569083 


4.569918 


4.569102 


5.479161 


5.479049 


30 


4.673213 


4.673208 


4.673227 


5.585293 


5.585280 


50 


4.703540 


4.703512 


4.703568 


5.623057 


5.622830 


70 




4.732207 


4.732279 


5.654646 


5.654747 


90 




4.775002 


4.775024 


5.695435 


5.695564 


110 




4.780718 


4.780809 


5.698338 


5.698447 


130 




4.786383 


4.786437 


5.704859 


5.704947 


150 




4.775124 


4.775231 


5.693268 


5.693447 
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The matrices initiating the sequence are obtained from 
!F° and T a by appropriately compressing their spectrum 
into the domain of convergence |l9| using 

x o = and x* = ^ t (8) 

J~ max min J~ min *' max 

where T m in and T max are approximate upper and lower 
bounds to the eigenvalues of J 70 . 

Recursion of the perturbed projection sequence is 
stopped when the change \Xf +1 — X°"\ becomes small. 
Having solved for the next Fock matrix deriva- 

tive F% + i is built, and the iteration continues until 
self-consistency, when the density matrix derivative 
and the desired property (e.g. the polarizability a ao — 
—2Tr[D^fib]) have reached a sufficient level of accuracy. 

We have implemented these methods in the Mon- 
doSCF suite of linear scaling quantum chemistry pro- 
grams |4^ . and performed polarizability calculations on 
a series of water clusters up to (H 2 0)150. These clus- 
ters were obtained by carving a spherical region out of 
a snapshot from a periodic classical molecular dynam- 
ics simulation of water at standard liquid density, and 
have been used previously in a number of scaling tests 

Hamuli. 

Calculations have been carried out at both the RHF/6- 
31G and RHF/6-31G** levels of theory and with both 
the GOOD and TIGHT thresholding parameter sets that 
control precision of the linear scaling algorithms, corre- 
sponding to matrix thresholds of 10~ 5 and 10 -6 , respec- 
tively. These calculations were carried out on a single 



Intel Xeon 2.4GHz processor running RedHat Linux 8.0 
and executables compiled with Portland Group Fortran 
Compiler pgf90 4.0-2 0. In Fig.HJ the total CPU time 
for the fifth CPSCF cycle (including build time for T a ; 
iterative construction of T> a and all intermediate steps 
including congruence transformation) is shown for the 
RHF/6-31G and RHF/6-31G** series of water clusters. 
Convergence of the CPSCF equations for these systems 
are typically achieved in about 10 cycles, independent 
of cluster size, basis set or matrix threshold. In Table 
[Q the corresponding average water cluster polarizabili- 
ties computed with the MondoSCF algorithms are listed 
and compared to the those obtained with the GAMESS 
quantum chemistry package at the RHF/6-31G level 
of theory. Figure [5] shows the magnitude of density ma- 
trix derivative atom-atom blocks as a function of atom- 
atom distance under global perturbation by a static elec- 
tric field. 

These results demonstrate an onset of linear scaling 
as early as 70 water molecules for properties with 4 dig- 
its of precision (RHF/6-31G at GOOD). While the CP- 
SCF equations must be solved iteratively with perturbed 
projection, the number of CPSCF cycles is ~ 10 when 
using DDIIS acceleration on well behaved systems. Us- 
ing an incomplete sparse linear algebra with thresholding 
is a particular advantage of the present implementation, 
as the small gap limit correctly leads to 0(N 3 ) algo- 
rithms while preserving accuracy. In contrast, methods 
that employ radial cutoffs incorrectly retain A^-scaling in 
this limit at the sacrifice of accuracy. However, the iter- 
ative approach proposed here for solution of the CPSCF 
equations is prone to the same instabilities encountered 
by the SCF equations in the small gap limit. 

We have presented a simple and efficient algorithm 
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FIG. 2: Magnitudes of atom-atom blocks of the RHF/6-31G 
density matrix derivative in the z direction with the sepa- 
ration of basis function centers for (-H2O)i50- The density 
matrix derivative has been converged using a TIGHT accuracy 
level (e.g. a matrix threshold of 10~ 6 au). 



numerical computations have been performed on com- 
puting resources located at this facility. 
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for solution of the coupled-perturbed self-consistent-field 
(CPSCF) equations in the context of spectral projection 
and the static electric polarizability. Unique features of 
the perturbed projection algorithm include linear scaling, 
simplicity, numerical stability and quadratic convergence 
in computation of the derivative density matrix. 

We have shown that the density matrix response is lo- 
cal upon global electric perturbation, corresponding to 
an approximate exponential decay of matrix elements. 
A similar exponential decay in the first order response 
corresponding to the local nuclear displacement has pre- 
viously been demonstrated by Ochsenfeld and Head- 
Gordon |2l| . These key observations are expected to hold 
generally for both local and global perturbations to insu- 
lating systems. The implication of these results are that 
the perturbed projection algorithm described in steps 
Ij5al5c|l and Eqs. (|6I7|I can be easily extended to the linear 
scaling computation of higher order response functions, 
DF and HF/DF models and a large class of static molec- 
ular properties such as the nuclear magnetic shielding 
tensor (NMR shift), indirect spin-spin coupling and the 
electronic g-tensor. We note also that the method is not 
unique to the TC2 generator or MondoSCF iV-scaling 
algorithms, but can be straightforwardly extended to 
other purification schemes such as TRS4 20] as well as 
other electronic structure programs. 
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